
% Create lookup file for checking which reactions and metabolites are
% present in which AGORA2 taxa
%
% OUTPUT
% lookupFilePath        Path to the generated lookup file
%
% Authors:  Tim Hensen, 2024

% Obtain paths to all species
panPath = 'C:\Users\mspg\Documents\Code\panSpecies';
% Step 1: Find paths to all agora2 taxa
panDir = what(panPath);

if isempty(panDir)
    error('panPath has not been found')
end

% Find paths to AGORA2 taxa
taxaPaths = string(strcat(panDir.path,filesep, panDir.mat));
taxa = erase(panDir.mat,'.mat');

% Remove duplicates
[taxa,idx] = unique(taxa);
taxaPaths = taxaPaths(idx);

tic
% Preallocate cell arrays to store strain data
taxaData = cell(1, length(taxa));

disp('Obtain reaction and metabolite content for each strain')
% Start parallel pool on all but two cores
if feature('numCores') > 5
    poolobj = gcp('nocreate');
    if isempty(poolobj)
        parpool(feature('numCores')-2)
    end
end
environment = getEnvironment();

disp('Obtain reaction and metabolite content for each strain')
parfor (i = 1:length(taxaPaths))
    % Restore environment
    restoreEnvironment(environment);
    changeCobraSolver('ibm_cplex', 'LP', 0, -1);

    disp(i)
    % Step 2: load each strain
    model = load(taxaPaths(i));
    model = model.(string(fieldnames(model)));

    % get all exchange reactions
    rxnIndex = find(contains(model.rxns,'EX_'));
    
    % Perform FVA
    [~,maxFlux] = guidedSim(model, model.rxns(rxnIndex));
    
    % Obtain produced metabolites
    rxnIndex(maxFlux<1e-6)=[];

    % Step 3: Create structure with strain as fieldname and reactions as field content
    taxonData{i}.taxon = taxa(i);
    taxonData{i}.rxns = model.rxns(rxnIndex);
end

% Concatenate strain data into database struct
database = struct;
for i = 1:length(taxa)
    fname = strcat('taxon_', string(i));
    database.(fname) = taxonData{i};
end

toc

% Step 4: Summarise structure to COBRA-like structure
disp('Process strain contents into an efficient data structure')

% Create structure
taxonMetProdCapacity = struct;

% Find all unique rxns, mets, and genes, and taxa
disp('Find the union of reactions, metabolites, and taxa')
taxonMetProdCapacity.rxns = obtainUnion(database, "rxns");
taxonMetProdCapacity.taxa = obtainUnion(database, "taxon");

% for each strain, indicate if a rxns, met, or gene is present
disp('Generate sparse matrix indicating which taxa contains which reactions and metabolites.')
taxonMetProdCapacity.rxnPresence = presenceMatrix(database, taxonMetProdCapacity,"rxns");
%
% Save taxonMetProdCapacity struct
disp('Save file in current directory')
lookupFilePath = [pwd filesep 'AGORA2_speciesMetCapacity.mat'];
save(lookupFilePath,'taxonMetProdCapacity')


function union_values = obtainUnion(database, type)
% Get fieldnames
field_names = fieldnames(database);
% Extract vectors from each field using cellfun
all_vectors_cell = cellfun(@(field) database.(field).(type), field_names, 'UniformOutput', false);
% Combine all vectors into one array
combined_vectors = cat(1, all_vectors_cell{:});
% Find the unique values in the combined array
union_values = unique(combined_vectors);
end


function matrix = presenceMatrix(database, taxonMetProdCapacity,type)
% Create matrix 
matrix = zeros(length(taxonMetProdCapacity.taxa), length(taxonMetProdCapacity.(type)));
% Check which reactions, metabolites, or genes in AGORA2 are present in the
% current strain
fname = fieldnames(database);
for i=1:length(fname)
    matrix(i,:) = matches(taxonMetProdCapacity.(type),database.(fname{i}).(type));
end
matrix = sparse(matrix);
end

